// QtCore
#include <QtCore/QDir>
#include <QtCore/QFile>
#include <QtCore/QTextCodec>
#include <QtCore/QDataStream>
#include <QtCore/QStringBuilder>
#include <QtCore/QCoreApplication>
// Локальные
#include "settings.h"
#include "math/math.h"
#include "linal/los.h"
#include "linal/bcg.h"
#include "solver_tpwft.h"
#include "fem/trianglefem.h"
#include "geometry/2d/vector.h"

//================================================================================================================================
//================================================================================================================================
SolverTPWFt::SolverTPWFt() : m_NetSearcher(NULL)
{
    //****************************************************************************************************************************
    // Создание и очистка директории для вывода результатов.
    //****************************************************************************************************************************
    m_OutputPath = Settings::instance()->setting("OutputPath").toString();
    if(!QDir::isAbsolutePath(m_OutputPath))
        m_OutputPath = QDir::cleanPath(m_OutputPath).prepend(qApp->applicationDirPath() + "/");

    if(!QDir(m_OutputPath).exists())
        QDir(qApp->applicationDirPath()).mkpath(m_OutputPath);
    if(!QDir(m_OutputPath + "/txt").exists())
        QDir(qApp->applicationDirPath()).mkpath(m_OutputPath + "/txt");
    // TODO: Очистка директории.

    //****************************************************************************************************************************
    // Считывание параметров временной сетки.
    //****************************************************************************************************************************
    m_NumOfTimeLayers = Settings::instance()->setting("Time/TimeLayers").toInt();
    m_TimeStep = Settings::instance()->setting("Time/TimeStep").toDouble();

    //****************************************************************************************************************************
    // Считывание физических и геометрических параметров задачи. Перевод скорости из см/ч в м/с.
    //****************************************************************************************************************************
    m_TCrystal = Settings::instance()->setting("Physics/CrystallizationTemperature").toDouble();
    m_TMax = Settings::instance()->setting("Physics/MaxTemperature").toDouble();
    m_TRod = Settings::instance()->setting("Physics/RodTemperature").toDouble();

    m_TGrad = Settings::instance()->setting("Physics/TemperatureGradient").toDouble();
    m_TGrad *= 100.;    //К/см -> К/м

    m_VRot = Settings::instance()->setting("Physics/AngularSpeed").toDouble();

    m_V = Settings::instance()->setting("Physics/CrucibleSpeed").toDouble();
    m_V /= 100.;        //см/ч -> м/ч
    m_V /= 3600.;       //м/ч -> м/с

    m_IsReversed = Settings::instance()->setting("Physics/IsReversed").toBool();
    m_ReverseStep = Settings::instance()->setting("Physics/ReverseStep").toInt();

    //****************************************************************************************************************************
    // Считывание параметров решателей.
    //****************************************************************************************************************************
    m_MaxiterSLAE = Settings::instance()->setting("Math/SlaeIters").toInt();
    m_MaxiterNonlin = Settings::instance()->setting("Math/NonlinIters").toInt();

    m_EpsSLAE = Settings::instance()->setting("Math/SlaeResidual").toDouble();
    m_EpsNonlin = Settings::instance()->setting("Math/NonlinResidual").toDouble();

    //****************************************************************************************************************************
    // Инициализация решателей.
    //****************************************************************************************************************************
    m_Solvers.append(new linal::LOS(m_MaxiterSLAE, m_EpsSLAE));
    m_Solvers.append(new linal::BCG(m_MaxiterSLAE, m_EpsSLAE));

    //****************************************************************************************************************************
    //****************************************************************************************************************************
    m_FrontHalfWidth = Settings::instance()->setting("Math/FrontHalfwidth").toDouble();
    m_MinCrystalHeight = Settings::instance()->setting("Math/MinCrystalHeight").toDouble();

    //****************************************************************************************************************************
    //****************************************************************************************************************************
    m_DiffStep = Settings::instance()->setting("Math/DiffStep").toDouble();

    //****************************************************************************************************************************
    //****************************************************************************************************************************
    m_W = Settings::instance()->setting("Math/RelaxFactor").toDouble();

    //****************************************************************************************************************************
    // Считывание сервисных настроек программы.
    //****************************************************************************************************************************
    m_OutputStep = Settings::instance()->setting("Solver/OtputStep").toInt();
}

//================================================================================================================================
//================================================================================================================================
bool SolverTPWFt::solveSLAE(const linal::CSRMatrix &A, const linal::Vector &f, linal::Vector &x)
{
    linal::Vector x_old = x;

    foreach(linal::SLAESolver* solver, m_Solvers)
    {
        x = x_old;
        if(solver->solveLU(A,f,x) < m_EpsSLAE)
            return true;
    }

    foreach(linal::SLAESolver* solver, m_Solvers)
    {
        x = x_old;
        if(solver->solveDi(A,f,x) < m_EpsSLAE)
            return true;
    }

    foreach(linal::SLAESolver* solver, m_Solvers)
    {
        x = x_old;
        if(solver->solve(A,f,x) < m_EpsSLAE)
            return true;
    }

    qWarning("SLAE wasn't solved!");
    return false;
}

//================================================================================================================================
//================================================================================================================================
void SolverTPWFt::solve()
{
    QTextStream cout(stdout);
    cout.setCodec(QTextCodec::codecForName("IBM866"));

    //****************************************************************************************************************************
    //****************************************************************************************************************************
    int resume_from = Settings::instance()->setting("Solver/ResumeFrom").toInt();

    if(resume_from != -1)
        readResumeInfo(resume_from);

    int first_layer = (resume_from == -1 ? 1 : resume_from + 1);
    for(int cur_layer = first_layer; cur_layer < m_NumOfTimeLayers; cur_layer++)
    {
        int nolin_iter(0);
        double nolin_eps(0.);

        do
        {
            formV();

            //Решаем уравнение для T
            m_A.zeroize();
            m_F.zeroize();
            formTSLAE(cur_layer);
            solveSLAE(m_A, m_F, m_T);
            m_T = (1 - m_W)*m_TOldNolin + m_W*m_T;

            recreateNet();
#ifndef QT_NO_DEBUG
//            writeResult(cur_layer);
#endif

            m_A.zeroize();
            m_F.zeroize();
            formVfSLAE(cur_layer);
            solveSLAE(m_A, m_F, m_Vf);
            m_Vf = (1 - m_W)*m_VfOldNolin + m_W*m_Vf;

            //Решаем уравнение для Omega
            m_A.zeroize();
            m_F.zeroize();
            formOmegaSLAE(cur_layer);
            solveSLAE(m_A, m_F, m_Omega);
            m_Omega = (1 - m_W)*m_OmegaOldNolin + m_W*m_Omega;

            //Решаем уравнение для Psi
            m_A.zeroize();
            m_F.zeroize();
            formPsiSLAE(cur_layer);
            solveSLAE(m_A, m_F, m_Psi);
            m_Psi = (1 - m_W)*m_PsiOldNolin + m_W*m_Psi;

            //Рассчитываем невязку по нелинейности
            double res_t = linal::Vector::normSquared(m_T - m_TOldNolin)/linal::Vector::normSquared(m_T);
            double res_vf = linal::Vector::normSquared(m_Vf - m_VfOldNolin)/linal::Vector::normSquared(m_Vf);
            double res_psi = linal::Vector::normSquared(m_Psi - m_PsiOldNolin)/linal::Vector::normSquared(m_Psi);
            double res_omega = linal::Vector::normSquared(m_Omega - m_OmegaOldNolin)/linal::Vector::normSquared(m_Omega);
            nolin_eps = qMax(qSqrt(res_t)/m_W, qMax(qSqrt(res_vf)/m_W, qMax(qSqrt(res_psi)/m_W, qSqrt(res_omega)/m_W)));

            m_TOldNolin = m_T;
            m_VfOldNolin = m_Vf;
            m_PsiOldNolin = m_Psi;
            m_OmegaOldNolin = m_Omega;
        }
        while(++nolin_iter < m_MaxiterNonlin && nolin_eps > m_EpsNonlin);

        m_TOldTime = m_T;
        m_VfOldTime = m_Vf;
        m_OmegaOldTime = m_Omega;

        if(cur_layer%m_OutputStep == 0)
            writeResult(cur_layer);
        cout << cur_layer << ": " << nolin_iter << " " << nolin_eps << endl;
    }
}

//================================================================================================================================
//================================================================================================================================
void SolverTPWFt::formV()
{
    m_Vr.resize(m_Net.numOfPoints());
    m_Vz.resize(m_Net.numOfPoints());

    //Находим r-компоненту скорости
    for(int i = 0; i < m_Vr.dimension(); i++)
    {
        if((!m_NetSearcher->pointBelongBoundary(i) && m_NetSearcher->pointMaterial(i) == 1) || m_NetSearcher->pointBelongBoundary(i,1,3))
        {
            geometry::_2d::Point p(m_Net.point(i).shifted(0, -m_DiffStep));
            int triangle = m_NetSearcher->findTriangle(m_Net, p);
            m_Vr[i] = - (fem::triangles::functionValue(triangle, m_Net, m_Psi, p) - m_Psi[i])/(-m_DiffStep);
        }
        else
            m_Vr[i] = 0;
    }

    //Находим z-компоненту скорости
    for(int i = 0; i < m_Vz.dimension(); i++)
    {
        if(!m_NetSearcher->pointBelongBoundary(i) && m_NetSearcher->pointMaterial(i) == 1)
        {
            geometry::_2d::Point p(m_Net.point(i).shifted(m_DiffStep, 0));
            int triangle = m_NetSearcher->findTriangle(m_Net, p);

            if(m_Net.point(i).x() == 0)
                m_Vz[i] = 2.*(fem::triangles::functionValue(triangle, m_Net, m_Psi, p) - m_Psi[i])/m_DiffStep;
            else
            {
                m_Vz[i] = (fem::triangles::functionValue(triangle, m_Net, m_Psi, p) - m_Psi[i])/m_DiffStep;
                m_Vz[i] += m_Psi[i]/m_Net.point(i).x();
            }
        }
        else
            m_Vz[i] = 0;
    }
}

//================================================================================================================================
//================================================================================================================================
void SolverTPWFt::tBC(int time_layer)
{
    for(int i = 0; i < m_TopTBC.size(); i++)
    {
        m_A.makeIdentity(m_TopTBC[i]);

        double time = time_layer*m_TimeStep;
        switch(m_TypeTBC[i])
        {
        case 0:
            m_F[m_TopTBC[i]] = m_TMax - m_V*m_TGrad*time;
            break;
        case 1:
            if(m_Net.point(m_TopTBC[i]).y() >= time*m_V)
                m_F[m_TopTBC[i]] = m_TMax;
            else
                m_F[m_TopTBC[i]] = m_TGrad*m_Net.point(m_TopTBC[i]).y() + m_TMax - m_V*m_TGrad*time;
            break;
        case 2:
            m_F[m_TopTBC[i]] = m_TRod;
            break;
        default:
            qFatal("Unknown BC type");
        }
    }
}

//================================================================================================================================
//================================================================================================================================
void SolverTPWFt::vfBC(int time_layer)
{
    Q_UNUSED(time_layer)

    for(int i = 0; i < m_TopVfBC.size(); i++)
    {
        m_A.makeIdentity(m_TopVfBC[i]);

        switch(m_TypeVfBC[i])
        {
        case 0:
            if(m_IsReversed == false || int(time_layer*m_TimeStep*m_VRot/(2*3.14))/m_ReverseStep%2 == 0)
                m_F[m_TopVfBC[i]] = m_VRot*m_Net.point(m_TopVfBC[i]).x();
            else
                m_F[m_TopVfBC[i]] = -m_VRot*m_Net.point(m_TopVfBC[i]).x();
            break;
        default:
            qFatal("Unknown BC type");
        }
    }
}

//================================================================================================================================
//================================================================================================================================
void SolverTPWFt::psiBC(int time_layer)
{
    Q_UNUSED(time_layer)

    for(int i = 0; i < m_TopPsiBC.size(); i++)
    {
        m_A.makeIdentity(m_TopPsiBC[i]);
        m_F[m_TopPsiBC[i]] = 0;
    }
}

//================================================================================================================================
//================================================================================================================================
void SolverTPWFt::omegaBC(int time_layer)
{
    Q_UNUSED(time_layer)

    for(int i = 0; i < m_TopOmegaBC.size(); i++)
    {
        m_A.makeIdentity(m_TopOmegaBC[i]);

        switch(m_TypeOmegaBC[i])
        {
        case 0: // Дно
        {
            double r = m_Net.point(m_TopOmegaBC[i]).x();
            double z = m_Net.point(m_TopOmegaBC[i]).y() + m_DiffStep;

            int triangle = m_NetSearcher->findTriangle(m_Net, geometry::_2d::Point(r, z));
            double psi = fem::triangles::functionValue(triangle, m_Net, m_Psi, geometry::_2d::Point(r, z));
            double omega = fem::triangles::functionValue(triangle, m_Net, m_OmegaOldNolin, geometry::_2d::Point(r, z));
            double woods = -omega/2 - 3*psi/(m_DiffStep*m_DiffStep);
            m_F[m_TopOmegaBC[i]] = woods;
        }
        break;
        case 1: // Ось симметрии
            m_F[m_TopOmegaBC[i]] = 0;
            break;
        case 2: // Боковая стенка
        {
            double r = m_Net.point(m_TopOmegaBC[i]).x() - m_DiffStep;
            double z = m_Net.point(m_TopOmegaBC[i]).y();

            int triangle = m_NetSearcher->findTriangle(m_Net, geometry::_2d::Point(r, z));
            double psi = fem::triangles::functionValue(triangle, m_Net, m_Psi, geometry::_2d::Point(r, z));
            double omega = fem::triangles::functionValue(triangle, m_Net, m_OmegaOldNolin, geometry::_2d::Point(r, z));
            double woods = -omega/2 - 3*psi/(m_DiffStep*m_DiffStep);
            m_F[m_TopOmegaBC[i]] = woods;
        }
        break;
        case 3: // Свободная поверхность
            m_F[m_TopOmegaBC[i]] = 0;
            break;
        case 4: // Фронт
        {
            int n = m_Front.indexOf(m_TopOmegaBC[i]);
            if(n ==0 || n == m_Front.size() - 1)
            {
                double r = m_Net.point(m_TopOmegaBC[i]).x();
                double z = m_Net.point(m_TopOmegaBC[i]).y() + m_DiffStep;

                int triangle = m_NetSearcher->findTriangle(m_Net, geometry::_2d::Point(r, z));
                double psi = fem::triangles::functionValue(triangle, m_Net, m_Psi, geometry::_2d::Point(r, z));
                double omega = fem::triangles::functionValue(triangle, m_Net, m_OmegaOldNolin, geometry::_2d::Point(r, z));
                double woods = -omega/2 - 3*psi/(m_DiffStep*m_DiffStep);
                m_F[m_TopOmegaBC[i]] = woods;
            }
            else
            {
                geometry::_2d::Vector direction = geometry::_2d::Vector::normal(geometry::_2d::Vector(m_Net.point(m_Front[n - 1]), m_Net.point(m_Front[n + 1])));
                direction.normalize();
                if(m_Net.point(m_Front[n]).shifted(direction.x(), direction.y()).y() - m_Net.point(m_Front[n]).y())
                    direction = -direction;

                double r = m_Net.point(m_TopOmegaBC[i]).x() + direction.x()*m_DiffStep;
                double z = m_Net.point(m_TopOmegaBC[i]).y() + direction.y()*m_DiffStep;

                int triangle = m_NetSearcher->findTriangle(m_Net, geometry::_2d::Point(r, z));
                double psi = fem::triangles::functionValue(triangle, m_Net, m_Psi, geometry::_2d::Point(r, z));
                double omega = fem::triangles::functionValue(triangle, m_Net, m_OmegaOldNolin, geometry::_2d::Point(r, z));
                double woods = -omega/2 - 3*psi/(m_DiffStep*m_DiffStep);
                m_F[m_TopOmegaBC[i]] = woods;
            }
        }
        case 5: // Наклонная стенка.
        {
            double r = m_Net.point(m_TopOmegaBC[i]).x() - m_DiffStep;
            //HACK: Потом подумать как сделать так, чтобы r заведомо не было равно нулю.
            if(r < 0) r = 0.;
            double z = m_Net.point(m_TopOmegaBC[i]).y() + m_DiffStep;

            int triangle = m_NetSearcher->findTriangle(m_Net, geometry::_2d::Point(r, z));
            double psi = fem::triangles::functionValue(triangle, m_Net, m_Psi, geometry::_2d::Point(r, z));
            double omega = fem::triangles::functionValue(triangle, m_Net, m_OmegaOldNolin, geometry::_2d::Point(r, z));
            double woods = -omega/2 - 3*psi/(2.*m_DiffStep*m_DiffStep);
            m_F[m_TopOmegaBC[i]] = woods;
        }
        break;
        default:
            qFatal("Unknown BC type");
        }
    }
}

//================================================================================================================================
//================================================================================================================================
void SolverTPWFt::formTSLAE(int time_layer)
{
    double density;
    double heat_capacity;
    double heat_conductivity;
    double D[3][3];             // Матрица, используемая для определения базисных функций
    double C[3][3][3];          // Локальная матрица массы
    double detD;                // Определитель матрицы D
    double abs_detD;            // Модуль определетиля матрицы D
    double ht;                  // Шаг по времени
    double r[3];                // r-координаты вершин КЭ
    double z[3];                // z-координаты вершин КЭ
    int top[3];                 // Индексы вершин КЭ

    ht = m_TimeStep;            // Задел на будущее, когда сетка может быть неравномерной.

    //****************************************************************************************************************************
    //Собираем конечноэлементную СЛАУ
    //****************************************************************************************************************************
    for(int cur_triangle = 0; cur_triangle < m_Net.numOfElements(); cur_triangle++)
    {
        density = m_Materials[m_Net.material(cur_triangle)].property(Material::Density);
        heat_capacity = m_Materials[m_Net.material(cur_triangle)].property(Material::HeatCapacity);
        heat_conductivity = m_Materials[m_Net.material(cur_triangle)].property(Material::HeatConductivity);

        for(int i = 0; i < 3; i++)
        {
            r[i] = m_Net.point(cur_triangle, i).x();
            z[i] = m_Net.point(cur_triangle, i).y();
            top[i] = m_Net.pointIndex(cur_triangle, i);
        }

        detD = (r[1] - r[0])*(z[2] - z[0]) - (r[2] - r[0])*(z[1] - z[0]);
        abs_detD = qAbs(detD);

        D[0][0] = (r[1]*z[2] - r[2]*z[1])/detD;
        D[0][1] = (z[1] - z[2])/detD;
        D[0][2] = (r[2] - r[1])/detD;
        D[1][0] = (r[2]*z[0] - r[0]*z[2])/detD;
        D[1][1] = (z[2] - z[0])/detD;
        D[1][2] = (r[0] - r[2])/detD;
        D[2][0] = (r[0]*z[1] - r[1]*z[0])/detD;
        D[2][1] = (z[0] - z[1])/detD;
        D[2][2] = (r[1] - r[0])/detD;

        C[0][0][0] = abs_detD/20.;      C[0][0][1] = abs_detD/60.;      C[0][0][2] = abs_detD/60.;
        C[0][1][0] = abs_detD/60.;      C[0][1][1] = abs_detD/60.;      C[0][1][2] = abs_detD/120.;
        C[0][2][0] = abs_detD/60.;      C[0][2][1] = abs_detD/120.;     C[0][2][2] = abs_detD/60.;

        C[1][0][0] = abs_detD/60.;      C[1][0][1] = abs_detD/60.;      C[1][0][2] = abs_detD/120.;
        C[1][1][0] = abs_detD/60.;      C[1][1][1] = abs_detD/20.;      C[1][1][2] = abs_detD/60.;
        C[1][2][0] = abs_detD/120.;     C[1][2][1] = abs_detD/60.;      C[1][2][2] = abs_detD/60.;

        C[2][0][0] = abs_detD/60.;      C[2][0][1] = abs_detD/120.;     C[2][0][2] = abs_detD/60.;
        C[2][1][0] = abs_detD/120.;     C[2][1][1] = abs_detD/60.;      C[2][1][2] = abs_detD/60.;
        C[2][2][0] = abs_detD/60.;      C[2][2][1] = abs_detD/60.;      C[2][2][2] = abs_detD/20.;

        for(int i = 0; i < 3; i++)
        {
            for(int j = 0; j < 3; j++)
            {
                double addition(0);
                addition += (r[0] + r[1] + r[2])*(D[i][1]*D[j][1] + D[i][2]*D[j][2])*heat_conductivity*abs_detD/6.;
                addition += (C[i][j][0]*r[0] + C[i][j][1]*r[1] + C[i][j][2]*r[2])*heat_capacity*density/ht;
                addition += (m_Vr[top[0]]*(r[0]*C[i][0][0] + r[1]*C[i][0][1] + r[2]*C[i][0][2]) +
                             m_Vr[top[1]]*(r[0]*C[i][1][0] + r[1]*C[i][1][1] + r[2]*C[i][1][2]) +
                             m_Vr[top[2]]*(r[0]*C[i][2][0] + r[1]*C[i][2][1] + r[2]*C[i][2][2]))*D[j][1]*heat_capacity*density;
                addition += (m_Vz[top[0]]*(r[0]*C[i][0][0] + r[1]*C[i][0][1] + r[2]*C[i][0][2]) +
                             m_Vz[top[1]]*(r[0]*C[i][1][0] + r[1]*C[i][1][1] + r[2]*C[i][1][2]) +
                             m_Vz[top[2]]*(r[0]*C[i][2][0] + r[1]*C[i][2][1] + r[2]*C[i][2][2]))*D[j][2]*heat_capacity*density;
                m_A(top[i], top[j]) += addition;
            }
        }

        for(int i = 0; i < 3; i++)
            for(int j = 0; j < 3; j++)
                for(int k = 0; k < 3; k++)
                    m_F[top[i]] += m_TOldTime[top[j]]*r[k]*C[i][j][k]*heat_capacity*density/ht;
    }

    //****************************************************************************************************************************
    //Учитываем краевые условия.
    //****************************************************************************************************************************
    tBC(time_layer);
}

//================================================================================================================================
//================================================================================================================================
void SolverTPWFt::formVfSLAE(int time_layer)
{
    double viscosity;
    double D[3][3];         //Матрица, используемая для определения базисных функций
    double C[3][3][3];      //Локальная матрица массы
    double detD;            //Определитель матрицы D
    double abs_detD;        //Модуль определетиля матрицы D
    double ht;              //Шаг по времени
    double r[3];            //r-координаты вершин КЭ
    double z[3];            //z-координаты вершин КЭ
    int top[3];             //Индексы вершин КЭ

    ht = m_TimeStep;        // Задел на будущее, когда сетка может быть неравномерной.

    //*******************************************************************************************************
    //Собираем конечноэлементную СЛАУ
    //*******************************************************************************************************
    for(int cur_triangle = 0; cur_triangle < m_Net.numOfElements(); cur_triangle++)
    {
        if(m_Materials[m_Net.material(cur_triangle)].state() == Material::Liquid)
        {
            viscosity = m_Materials[m_Net.material(cur_triangle)].property(Material::Viscosity);

            for(int i = 0; i < 3; i++)
            {
                r[i] = m_Net.point(cur_triangle, i).x();
                z[i] = m_Net.point(cur_triangle, i).y();
                top[i] = m_Net.pointIndex(cur_triangle, i);
            }

            detD = (r[1] - r[0])*(z[2] - z[0]) - (r[2] - r[0])*(z[1] - z[0]);
            abs_detD = qAbs(detD);

            D[0][0] = (r[1]*z[2] - r[2]*z[1])/detD;
            D[0][1] = (z[1] - z[2])/detD;
            D[0][2] = (r[2] - r[1])/detD;
            D[1][0] = (r[2]*z[0] - r[0]*z[2])/detD;
            D[1][1] = (z[2] - z[0])/detD;
            D[1][2] = (r[0] - r[2])/detD;
            D[2][0] = (r[0]*z[1] - r[1]*z[0])/detD;
            D[2][1] = (z[0] - z[1])/detD;
            D[2][2] = (r[1] - r[0])/detD;

            C[0][0][0] = abs_detD/20.;      C[0][0][1] = abs_detD/60.;      C[0][0][2] = abs_detD/60.;
            C[0][1][0] = abs_detD/60.;      C[0][1][1] = abs_detD/60.;      C[0][1][2] = abs_detD/120.;
            C[0][2][0] = abs_detD/60.;      C[0][2][1] = abs_detD/120.;     C[0][2][2] = abs_detD/60.;

            C[1][0][0] = abs_detD/60.;      C[1][0][1] = abs_detD/60.;      C[1][0][2] = abs_detD/120.;
            C[1][1][0] = abs_detD/60.;      C[1][1][1] = abs_detD/20.;      C[1][1][2] = abs_detD/60.;
            C[1][2][0] = abs_detD/120.;     C[1][2][1] = abs_detD/60.;      C[1][2][2] = abs_detD/60.;

            C[2][0][0] = abs_detD/60.;      C[2][0][1] = abs_detD/120.;     C[2][0][2] = abs_detD/60.;
            C[2][1][0] = abs_detD/120.;     C[2][1][1] = abs_detD/60.;      C[2][1][2] = abs_detD/60.;
            C[2][2][0] = abs_detD/60.;      C[2][2][1] = abs_detD/60.;      C[2][2][2] = abs_detD/20.;

            for(int i = 0; i < 3; i++)
            {
                for(int j = 0; j < 3; j++)
                {
                    double addition(0);
                    addition += (r[0] + r[1] + r[2])*(D[i][1]*D[j][1] + D[i][2]*D[j][2])*viscosity*abs_detD/6.;
                    addition += (C[i][j][0]*r[0] + C[i][j][1]*r[1] + C[i][j][2]*r[2])/ht;
                    addition += (m_Vr[top[0]]*(r[0]*C[i][0][0] + r[1]*C[i][0][1] + r[2]*C[i][0][2]) +
                                 m_Vr[top[1]]*(r[0]*C[i][1][0] + r[1]*C[i][1][1] + r[2]*C[i][1][2]) +
                                 m_Vr[top[2]]*(r[0]*C[i][2][0] + r[1]*C[i][2][1] + r[2]*C[i][2][2]))*D[j][1];
                    addition += (m_Vz[top[0]]*(r[0]*C[i][0][0] + r[1]*C[i][0][1] + r[2]*C[i][0][2]) +
                                 m_Vz[top[1]]*(r[0]*C[i][1][0] + r[1]*C[i][1][1] + r[2]*C[i][1][2]) +
                                 m_Vz[top[2]]*(r[0]*C[i][2][0] + r[1]*C[i][2][1] + r[2]*C[i][2][2]))*D[j][2];
                    addition += m_Vr[top[0]]*C[i][j][0];
                    addition += m_Vr[top[1]]*C[i][j][1];
                    addition += m_Vr[top[2]]*C[i][j][2];

                    m_A(top[i], top[j]) += addition;
                }
            }

            for(int i = 0; i < 3; i++)
                for(int j = 0; j < 3; j++)
                    m_F[top[i]] += m_VfOldTime[top[j]]*(r[0]*C[i][j][0] + r[1]*C[i][j][1] + r[2]*C[i][j][2])/ht;


            for(int i = 0; i < 3; i++)
                for(int j = 0; j < 3; j++)
                    if(r[j] != 0)
                        m_F[top[i]] -= m_VfOldNolin[top[j]]*(C[i][j][0] + C[i][j][1] + C[i][j][2])*viscosity/r[j];
                    else
                    {
                        geometry::_2d::Point p(m_Net.point(top[j]).shifted(m_DiffStep, 0));
                        int triangle = m_NetSearcher->findTriangle(m_Net, p);
                        double vf_diff = (fem::triangles::functionValue(triangle, m_Net, m_VfOldNolin, p) - m_VfOldNolin[top[j]])/m_DiffStep;
                        m_F[top[i]] -= vf_diff*(C[i][j][0] + C[i][j][1] + C[i][j][2])*viscosity;
                    }
        }
    }

    //*******************************************************************************************************
    //Учитываем краевые условия
    //*******************************************************************************************************
    vfBC(time_layer);

    //*******************************************************************************************************
    //Заполняем матрицу для нерасплава
    //*******************************************************************************************************
    for(int i = 0; i < m_A.dimension(); ++i)
        if(qAbs(m_A(i,i)) < DBL_EPSILON)
        {
            m_A(i, i) = 1.;
            if(m_IsReversed == false || int(time_layer*m_TimeStep*m_VRot/(2*3.14))/m_ReverseStep%2 == 0)
                m_F[i] = m_VRot*m_Net.point(i).x();
            else
                m_F[i] = -m_VRot*m_Net.point(i).x();
        }
}

//================================================================================================================================
//================================================================================================================================
void SolverTPWFt::formPsiSLAE(int time_layer)
{
    double D[3][3];         // Матрица, используемая для определения базисных функций
    double C[3][3][3];      // Локальная матрица массы
    double detD;            // Определитель матрицы D
    double abs_detD;        // Модуль определетиля матрицы D
    double r[3];            // r-координаты вершин КЭ
    double z[3];            // z-координаты вершин КЭ
    int top[3];             // Индексы вершин КЭ

    //****************************************************************************************************************************
    //Собираем конечноэлементную СЛАУ
    //****************************************************************************************************************************
    for(int cur_triangle = 0; cur_triangle < m_Net.numOfElements(); cur_triangle++)
    {
        if(m_Materials[m_Net.material(cur_triangle)].state() != Material::Liquid)
        {
            m_A(m_Net.pointIndex(cur_triangle, 0), m_Net.pointIndex(cur_triangle, 0)) = 1.;
            m_A(m_Net.pointIndex(cur_triangle, 1), m_Net.pointIndex(cur_triangle, 1)) = 1.;
            m_A(m_Net.pointIndex(cur_triangle, 2), m_Net.pointIndex(cur_triangle, 2)) = 1.;
        }
        else
        {
            for(int i = 0; i < 3; i++)
            {
                r[i] = m_Net.point(cur_triangle, i).x();
                z[i] = m_Net.point(cur_triangle, i).y();
                top[i] = m_Net.pointIndex(cur_triangle, i);
            }

            detD = (r[1] - r[0])*(z[2] - z[0]) - (r[2] - r[0])*(z[1] - z[0]);
            abs_detD = qAbs(detD);

            D[0][0] = (r[1]*z[2] - r[2]*z[1])/detD;
            D[0][1] = (z[1] - z[2])/detD;
            D[0][2] = (r[2] - r[1])/detD;
            D[1][0] = (r[2]*z[0] - r[0]*z[2])/detD;
            D[1][1] = (z[2] - z[0])/detD;
            D[1][2] = (r[0] - r[2])/detD;
            D[2][0] = (r[0]*z[1] - r[1]*z[0])/detD;
            D[2][1] = (z[0] - z[1])/detD;
            D[2][2] = (r[1] - r[0])/detD;

            C[0][0][0] = abs_detD/20.;      C[0][0][1] = abs_detD/60.;      C[0][0][2] = abs_detD/60.;
            C[0][1][0] = abs_detD/60.;      C[0][1][1] = abs_detD/60.;      C[0][1][2] = abs_detD/120.;
            C[0][2][0] = abs_detD/60.;      C[0][2][1] = abs_detD/120.;     C[0][2][2] = abs_detD/60.;

            C[1][0][0] = abs_detD/60.;      C[1][0][1] = abs_detD/60.;      C[1][0][2] = abs_detD/120.;
            C[1][1][0] = abs_detD/60.;      C[1][1][1] = abs_detD/20.;      C[1][1][2] = abs_detD/60.;
            C[1][2][0] = abs_detD/120.;     C[1][2][1] = abs_detD/60.;      C[1][2][2] = abs_detD/60.;

            C[2][0][0] = abs_detD/60.;      C[2][0][1] = abs_detD/120.;     C[2][0][2] = abs_detD/60.;
            C[2][1][0] = abs_detD/120.;     C[2][1][1] = abs_detD/60.;      C[2][1][2] = abs_detD/60.;
            C[2][2][0] = abs_detD/60.;      C[2][2][1] = abs_detD/60.;      C[2][2][2] = abs_detD/20.;

            for(int i = 0; i < 3; i++)
                for(int j = 0; j < 3; j++)
                    m_A(top[i], top[j]) += (r[0] + r[1] + r[2])*(D[i][1]*D[j][1] + D[i][2]*D[j][2])*abs_detD/6.;

            for(int i = 0; i < 3; i++)
                for(int j = 0; j < 3; j++)
                    m_F[top[i]] += m_Omega[top[j]]*(r[0]*C[i][j][0] + r[1]*C[i][j][1] + r[2]*C[i][j][2]);

            for(int i = 0; i < 3; i++)
                for(int j = 0; j < 3; j++)
                    if(r[j] != 0)
                        m_F[top[i]] -= m_PsiOldNolin[top[j]]*(C[i][j][0] + C[i][j][1] + C[i][j][2])/r[j];
                    else
                    {
                        geometry::_2d::Point p(m_Net.point(top[j]).shifted(m_DiffStep, 0));
                        int triangle = m_NetSearcher->findTriangle(m_Net, p);
                        double psi_diff = (fem::triangles::functionValue(triangle, m_Net, m_PsiOldNolin, p) - m_PsiOldNolin[top[j]])/m_DiffStep;

                        m_F[top[i]] -= psi_diff*(C[i][j][0] + C[i][j][1] + C[i][j][2]);
                    }
        }
    }

    //*******************************************************************************************************
    //Учитываем краевые условия
    //*******************************************************************************************************
    psiBC(time_layer);
}

//================================================================================================================================
//================================================================================================================================
void SolverTPWFt::formOmegaSLAE(int time_layer)
{
    double viscosity;
    double expansiveness;
    double D[3][3];         // Матрица, используемая для определения базисных функций
    double C[3][3][3];      // Локальная матрица массы
    double detD;            // Определитель матрицы D
    double abs_detD;        // Модуль определетиля матрицы D
    double ht;              // Шаг по времени
    double r[3];            // r-координаты вершин КЭ
    double z[3];            // z-координаты вершин КЭ
    int top[3];             // Индексы вершин КЭ

    ht = m_TimeStep;        // Задел на будущее, когда сетка может быть неравномерной.

    //*******************************************************************************************************
    //Собираем конечноэлементную СЛАУ
    //*******************************************************************************************************
    for(int cur_triangle = 0; cur_triangle < m_Net.numOfElements(); cur_triangle++)
    {
        if(m_Materials[m_Net.material(cur_triangle)].state() != Material::Liquid)
        {
            m_A(m_Net.pointIndex(cur_triangle, 0), m_Net.pointIndex(cur_triangle, 0)) = 1.;
            m_A(m_Net.pointIndex(cur_triangle, 1), m_Net.pointIndex(cur_triangle, 1)) = 1.;
            m_A(m_Net.pointIndex(cur_triangle, 2), m_Net.pointIndex(cur_triangle, 2)) = 1.;
        }
        else
        {
            viscosity = m_Materials[m_Net.material(cur_triangle)].property(Material::Viscosity);
            expansiveness = m_Materials[m_Net.material(cur_triangle)].property(Material::Expansiveness);

            for(int i = 0; i < 3; i++)
            {
                r[i] = m_Net.point(cur_triangle, i).x();
                z[i] = m_Net.point(cur_triangle, i).y();
                top[i] = m_Net.pointIndex(cur_triangle, i);
            }

            detD = (r[1] - r[0])*(z[2] - z[0]) - (r[2] - r[0])*(z[1] - z[0]);
            abs_detD = qAbs(detD);

            D[0][0] = (r[1]*z[2] - r[2]*z[1])/detD;
            D[0][1] = (z[1] - z[2])/detD;
            D[0][2] = (r[2] - r[1])/detD;
            D[1][0] = (r[2]*z[0] - r[0]*z[2])/detD;
            D[1][1] = (z[2] - z[0])/detD;
            D[1][2] = (r[0] - r[2])/detD;
            D[2][0] = (r[0]*z[1] - r[1]*z[0])/detD;
            D[2][1] = (z[0] - z[1])/detD;
            D[2][2] = (r[1] - r[0])/detD;

            C[0][0][0] = abs_detD/20.;      C[0][0][1] = abs_detD/60.;      C[0][0][2] = abs_detD/60.;
            C[0][1][0] = abs_detD/60.;      C[0][1][1] = abs_detD/60.;      C[0][1][2] = abs_detD/120.;
            C[0][2][0] = abs_detD/60.;      C[0][2][1] = abs_detD/120.;     C[0][2][2] = abs_detD/60.;

            C[1][0][0] = abs_detD/60.;      C[1][0][1] = abs_detD/60.;      C[1][0][2] = abs_detD/120.;
            C[1][1][0] = abs_detD/60.;      C[1][1][1] = abs_detD/20.;      C[1][1][2] = abs_detD/60.;
            C[1][2][0] = abs_detD/120.;     C[1][2][1] = abs_detD/60.;      C[1][2][2] = abs_detD/60.;

            C[2][0][0] = abs_detD/60.;      C[2][0][1] = abs_detD/120.;     C[2][0][2] = abs_detD/60.;
            C[2][1][0] = abs_detD/120.;     C[2][1][1] = abs_detD/60.;      C[2][1][2] = abs_detD/60.;
            C[2][2][0] = abs_detD/60.;      C[2][2][1] = abs_detD/60.;      C[2][2][2] = abs_detD/20.;

            for(int i = 0; i < 3; i++)
            {
                for(int j = 0; j < 3; j++)
                {
                    double addition(0);
                    addition += (r[0] + r[1] + r[2])*(D[i][1]*D[j][1] + D[i][2]*D[j][2])*viscosity*abs_detD/6.;
                    addition += (C[i][j][0]*r[0] + C[i][j][1]*r[1] + C[i][j][2]*r[2])/ht;
                    addition += (m_Vr[top[0]]*(r[0]*C[i][0][0] + r[1]*C[i][0][1] + r[2]*C[i][0][2]) +
                                 m_Vr[top[1]]*(r[0]*C[i][1][0] + r[1]*C[i][1][1] + r[2]*C[i][1][2]) +
                                 m_Vr[top[2]]*(r[0]*C[i][2][0] + r[1]*C[i][2][1] + r[2]*C[i][2][2]))*D[j][1];
                    addition += (m_Vz[top[0]]*(r[0]*C[i][0][0] + r[1]*C[i][0][1] + r[2]*C[i][0][2]) +
                                 m_Vz[top[1]]*(r[0]*C[i][1][0] + r[1]*C[i][1][1] + r[2]*C[i][1][2]) +
                                 m_Vz[top[2]]*(r[0]*C[i][2][0] + r[1]*C[i][2][1] + r[2]*C[i][2][2]))*D[j][2];
                    addition -= m_Vr[top[0]]*C[i][j][0];
                    addition -= m_Vr[top[1]]*C[i][j][1];
                    addition -= m_Vr[top[2]]*C[i][j][2];

                    m_A(top[i], top[j]) += addition;
                }
            }

            for(int i = 0; i < 3; i++)
                for(int j = 0; j < 3; j++)
                    for(int k = 0; k < 3; k++)
                        m_F[top[i]] -= m_T[top[j]]*r[k]*(C[i][k][0] + C[i][k][1] + C[i][k][2])*D[j][1]*expansiveness*9.8;

            for(int i = 0; i < 3; i++)
                for(int j = 0; j < 3; j++)
                    m_F[top[i]] += m_OmegaOldTime[top[j]]*(r[0]*C[i][j][0] + r[1]*C[i][j][1] + r[2]*C[i][j][2])/ht;

            for(int i = 0; i < 3; i++)
                for(int j = 0; j < 3; j++)
                    for(int k = 0; k < 3; k++)
                        m_F[top[i]] += 2*m_Vf[top[j]]*m_Vf[top[k]]*(C[i][j][0] + C[i][j][1] + C[i][j][2])*D[k][2];

            for(int i = 0; i < 3; i++)
                for(int j = 0; j < 3; j++)
                    if(r[j] != 0)
                        m_F[top[i]] -= m_OmegaOldNolin[top[j]]*(C[i][j][0] + C[i][j][1] + C[i][j][2])*viscosity/r[j];
                    else
                    {
                        geometry::_2d::Point p(m_Net.point(top[j]).shifted(m_DiffStep, 0));
                        int triangle = m_NetSearcher->findTriangle(m_Net, p);
                        double omega_diff = (fem::triangles::functionValue(triangle, m_Net, m_OmegaOldNolin, p) - m_OmegaOldNolin[top[j]])/(m_DiffStep);
                        m_F[top[i]] -= omega_diff*(C[i][j][0] + C[i][j][1] + C[i][j][2])*viscosity;
                    }
        }
    }

    //*******************************************************************************************************
    //Учитываем краевые условия
    //*******************************************************************************************************
    omegaBC(time_layer);
}

//================================================================================================================================
//Вывод результата в файл.
//================================================================================================================================
void SolverTPWFt::writeResult(int time_layer)
{
    QFile file_bin;
    QFile file_txt;
    QDataStream ostream_bin(&file_bin);
    QTextStream ostream_txt(&file_txt);
    ostream_bin.setVersion(QDataStream::Qt_4_7);
    ostream_txt.setRealNumberPrecision(12);

    // Количество десятичных разрядов в числе, равном количеству временных слоев.
    // Не все программы способны "правильно" упорядочивать файлы, в именах которых встречаются числа с разным количеством разрядов.
    // Так, например, одни программы считают, что файл 10.dat должен распологаться раньше файла 100.dat, в то время как другие
    // программы наоборот располагают файл 100.dat перед файлом 10.dat. Единственный способ заставить все программы одинаково
    // работать с подобными файлами - выравнять количество разрядов в числах, присутствующих в именах файлов, добавляя при
    // необходимости незначащие нули перед числом. Т.е. файл 010.dat любой программой будет расположен перед файлом 100.dat.
    // В данном приложении, в имени файлов указываются номер временного слоя и, сл-но, для того, чтобы знать нужно ли добавлять
    // незначащие нули, и если да, то сколько, необходимо знать число разрядов в максимально возможном номере временного слоя.
    int num_of_orders = QString::number(m_NumOfTimeLayers).count();

    //****************************************************************************************************************************
    //Выводим конечноэлементую сетку.
    //****************************************************************************************************************************
    QString filename_nvtr = QString("/%1nvtr.dat").arg(time_layer, num_of_orders, 10, QChar('0')).prepend(m_OutputPath);
    QString filename_rz = QString("/%1rz.dat").arg(time_layer, num_of_orders, 10, QChar('0')).prepend(m_OutputPath);
    m_Net.writeToFiles(filename_rz, filename_nvtr);

    file_txt.setFileName(QString("/txt/%1nvtr.dat").arg(time_layer, num_of_orders, 10, QChar('0')).prepend(m_OutputPath));
    file_txt.open(QFile::WriteOnly | QFile::Text);
    for(int i = 0; i < m_Net.numOfElements(); ++i)
    {
        for(int j = 0; j < 3; ++j)
            ostream_txt << m_Net.pointIndex(i, j) << " ";
        ostream_txt << m_Net.material(i) << endl;
    }
    file_txt.close();

    file_txt.setFileName(QString("/txt/%1rz.dat").arg(time_layer, num_of_orders, 10, QChar('0')).prepend(m_OutputPath));
    file_txt.open(QFile::WriteOnly | QFile::Text);
    for(int i = 0; i < m_Net.numOfPoints(); ++i)
        ostream_txt << m_Net.point(i) << endl;
    file_txt.close();

    //****************************************************************************************************************************
    //Выводим значения температуры.
    //****************************************************************************************************************************
    file_bin.setFileName(QString("/%1T.dat").arg(time_layer, num_of_orders, 10, QChar('0')).prepend(m_OutputPath));
    file_bin.open(QFile::WriteOnly);
    for(int i = 0; i < m_Net.numOfPoints(); i++)
        ostream_bin << m_Net.point(i).x() << m_Net.point(i).y() << m_T[i];
    file_bin.close();

    file_txt.setFileName(QString("/txt/%1T.dat").arg(time_layer, num_of_orders, 10, QChar('0')).prepend(m_OutputPath));
    file_txt.open(QFile::WriteOnly | QFile::Text);
    for(int i = 0; i < m_Net.numOfPoints(); i++)
        ostream_txt << m_Net.point(i) << " " << m_T[i] << endl;
    file_txt.close();

    //************************************************************************************************************************
    //Выводим значения функции тока.
    //************************************************************************************************************************
    file_bin.setFileName(QString("/%1Psi.dat").arg(time_layer, num_of_orders, 10, QChar('0')).prepend(m_OutputPath));
    file_bin.open(QFile::WriteOnly);
    for(int i = 0; i < m_Net.numOfPoints(); i++)
        ostream_bin << m_Net.point(i).x() << m_Net.point(i).y() << m_Psi[i];
    file_bin.close();

    file_txt.setFileName(QString("/txt/%1Psi.dat").arg(time_layer, num_of_orders, 10, QChar('0')).prepend(m_OutputPath));
    file_txt.open(QFile::WriteOnly | QFile::Text);
    for(int i = 0; i < m_Net.numOfPoints(); i++)
        ostream_txt << m_Net.point(i) << " " << m_Psi[i] << endl;
    file_txt.close();

    //************************************************************************************************************************
    //Выводим значения вихря.
    //************************************************************************************************************************
    file_bin.setFileName(QString("/%1Omega.dat").arg(time_layer, num_of_orders, 10, QChar('0')).prepend(m_OutputPath));
    file_bin.open(QFile::WriteOnly);
    for(int i = 0; i < m_Net.numOfPoints(); i++)
        ostream_bin << m_Net.point(i).x() << m_Net.point(i).y() << m_Omega[i];
    file_bin.close();

    file_txt.setFileName(QString("/txt/%1Omega.dat").arg(time_layer, num_of_orders, 10, QChar('0')).prepend(m_OutputPath));
    file_txt.open(QFile::WriteOnly | QFile::Text);
    for(int i = 0; i < m_Net.numOfPoints(); i++)
        ostream_txt << m_Net.point(i) << " " << m_Omega[i] << endl;
    file_txt.close();

    //************************************************************************************************************************
    //Выводим значения r-компоненты скорости.
    //************************************************************************************************************************
    file_bin.setFileName(QString("/%1Vr.dat").arg(time_layer, num_of_orders, 10, QChar('0')).prepend(m_OutputPath));
    file_bin.open(QFile::WriteOnly);
    for(int i = 0; i < m_Net.numOfPoints(); i++)
        ostream_bin << m_Net.point(i).x() << m_Net.point(i).y() << m_Vr[i];
    file_bin.close();

    file_txt.setFileName(QString("/txt/%1Vr.dat").arg(time_layer, num_of_orders, 10, QChar('0')).prepend(m_OutputPath));
    file_txt.open(QFile::WriteOnly | QFile::Text);
    for(int i = 0; i < m_Net.numOfPoints(); i++)
        ostream_txt << m_Net.point(i) << " " << m_Vr[i] << endl;
    file_txt.close();

    //************************************************************************************************************************
    //выводим значения z-компоненты скорости.
    //************************************************************************************************************************
    file_bin.setFileName(QString("/%1Vz.dat").arg(time_layer, num_of_orders, 10, QChar('0')).prepend(m_OutputPath));
    file_bin.open(QFile::WriteOnly);
    for(int i = 0; i < m_Net.numOfPoints(); i++)
        ostream_bin << m_Net.point(i).x() << m_Net.point(i).y() << m_Vz[i];
    file_bin.close();

    file_txt.setFileName(QString("/txt/%1Vz.dat").arg(time_layer, num_of_orders, 10, QChar('0')).prepend(m_OutputPath));
    file_txt.open(QFile::WriteOnly | QFile::Text);
    for(int i = 0; i < m_Net.numOfPoints(); i++)
        ostream_txt << m_Net.point(i) << " " << m_Vz[i] << endl;
    file_txt.close();

    //************************************************************************************************************************
    //выводим значения f-компоненты скорости.
    //************************************************************************************************************************
    file_bin.setFileName(QString("/%1Vf.dat").arg(time_layer, num_of_orders, 10, QChar('0')).prepend(m_OutputPath));
    file_bin.open(QFile::WriteOnly);
    for(int i = 0; i < m_Net.numOfPoints(); i++)
        ostream_bin << m_Net.point(i).x() << m_Net.point(i).y() << m_Vf[i];
    file_bin.close();

    file_txt.setFileName(QString("/txt/%1Vf.dat").arg(time_layer, num_of_orders, 10, QChar('0')).prepend(m_OutputPath));
    file_txt.open(QFile::WriteOnly | QFile::Text);
    for(int i = 0; i < m_Net.numOfPoints(); i++)
        ostream_txt << m_Net.point(i) << " " << m_Vf[i] << endl;
    file_txt.close();

    //************************************************************************************************************************
    //выводим информацию для восстановления решения.
    //************************************************************************************************************************
    file_bin.setFileName(QString("/%1bc.dat").arg(time_layer, num_of_orders, 10, QChar('0')).prepend(m_OutputPath));
    file_bin.open(QIODevice::WriteOnly);
    ostream_bin << m_TopTBC.size();
    for(int i = 0; i < m_TopTBC.size(); i++)
        ostream_bin << m_TopTBC[i] << m_TypeTBC[i];

    ostream_bin << m_TopVfBC.size();
    for(int i = 0; i < m_TopVfBC.size(); i++)
        ostream_bin << m_TopVfBC[i] << m_TypeVfBC[i];

    ostream_bin << m_TopPsiBC.size();
    for(int i = 0; i < m_TopPsiBC.size(); i++)
        ostream_bin << m_TopPsiBC[i] << m_TypePsiBC[i];

    ostream_bin << m_TopOmegaBC.size();
    for(int i = 0; i < m_TopOmegaBC.size(); i++)
        ostream_bin << m_TopOmegaBC[i] << m_TypeOmegaBC[i];
    file_bin.close();

    file_bin.setFileName(QString("/%1state.dat").arg(time_layer, num_of_orders, 10, QChar('0')).prepend(m_OutputPath));
    file_bin.open(QIODevice::WriteOnly);
    ostream_bin << m_State;
    file_bin.close();
}

//================================================================================================================================
//Инициализация структур данных, необходимая для продолжения расчетов с указанного места.
//================================================================================================================================
void SolverTPWFt::readResumeInfo(int time_layer)
{
    // Количество десятичных разрядов в числе, равном количеству временных слоев.
    // Не все программы способны "правильно" упорядочивать файлы, в именах которых встречаются числа с разным количеством разрядов.
    // Так, например, одни программы считают, что файл 10.dat должен распологаться раньше файла 100.dat, в то время как другие
    // программы наоборот располагают файл 100.dat перед файлом 10.dat. Единственный способ заставить все программы одинаково
    // работать с подобными файлами - выравнять количество разрядов в числах, присутствующих в именах файлов, добавляя при
    // необходимости незначащие нули перед числом. Т.е. файл 010.dat любой программой будет расположен перед файлом 100.dat.
    // В данном приложении, в имени файлов указываются номер временного слоя и, сл-но, для того, чтобы знать нужно ли добавлять
    // незначащие нули, и если да, то сколько, необходимо знать число разрядов в максимально возможном номере временного слоя.
    int num_of_orders = QString::number(m_NumOfTimeLayers).count();

    QString filename_rz = QString("/%1rz.dat").arg(time_layer, num_of_orders, 10, QChar('0')).prepend(m_OutputPath);
    QString filename_nvtr = QString("/%1nvtr.dat").arg(time_layer, num_of_orders, 10, QChar('0')).prepend(m_OutputPath);

    QString filename_t = QString("/%1T.dat").arg(time_layer, num_of_orders, 10, QChar('0')).prepend(m_OutputPath);
    QString filename_vf = QString("/%1Vf.dat").arg(time_layer, num_of_orders, 10, QChar('0')).prepend(m_OutputPath);
    QString filename_psi = QString("/%1Psi.dat").arg(time_layer, num_of_orders, 10, QChar('0')).prepend(m_OutputPath);
    QString filename_omega = QString("/%1Omega.dat").arg(time_layer, num_of_orders, 10, QChar('0')).prepend(m_OutputPath);

    QString filename_bc = QString("/%1bc.dat").arg(time_layer, num_of_orders, 10, QChar('0')).prepend(m_OutputPath);
    QString filename_state = QString("/%1state.dat").arg(time_layer, num_of_orders, 10, QChar('0')).prepend(m_OutputPath);

    //****************************************************************************************************************************
    //****************************************************************************************************************************
    m_Net.readFromFiles(filename_rz, filename_nvtr);
    m_NetSearcher->analyzeNet(m_Net);

    //****************************************************************************************************************************
    //****************************************************************************************************************************
    int data_size;
    int bc_top, bc_type;

    QFile file(filename_bc);
    QDataStream istream(&file);
    istream.setVersion(QDataStream::Qt_4_7);
    file.open(QIODevice::ReadOnly);

    istream >> data_size;
    m_TopTBC.clear();
    m_TypeTBC.clear();
    for(int i = 0; i < data_size; i++)
    {
        istream >> bc_top >> bc_type;
        m_TopTBC.append(bc_top);
        m_TypeTBC.append(bc_type);
    }

    istream >> data_size;
    m_TopVfBC.clear();
    m_TypeVfBC.clear();
    for(int i = 0; i < data_size; i++)
    {
        istream >> bc_top >> bc_type;
        m_TopVfBC.append(bc_top);
        m_TypeVfBC.append(bc_type);
    }

    istream >> data_size;
    m_TopPsiBC.clear();
    m_TypePsiBC.clear();
    for(int i = 0; i < data_size; i++)
    {
        istream >> bc_top >> bc_type;
        m_TopPsiBC.append(bc_top);
        m_TypePsiBC.append(bc_type);
    }

    istream >> data_size;
    m_TopOmegaBC.clear();
    m_TypeOmegaBC.clear();
    for(int i = 0; i < data_size; i++)
    {
        istream >> bc_top >> bc_type;
        m_TopOmegaBC.append(bc_top);
        m_TypeOmegaBC.append(bc_type);
    }

    file.close();

    //****************************************************************************************************************************
    //****************************************************************************************************************************
    int state;

    file.setFileName(filename_state);
    file.open(QIODevice::ReadOnly);
    istream >> state;
    switch(state)
    {
    case 0: m_State = SolverTPWFt::NotStarted; break;
    case 1: m_State = SolverTPWFt::JustStarted; break;
    case 2: m_State = SolverTPWFt::InProgress; break;
    case 3: m_State = SolverTPWFt::Finished; break;
    default: qFatal("Unknown crystallization state");
    }
    file.close();

    //****************************************************************************************************************************
    //****************************************************************************************************************************
    m_T.readFromFile(filename_t);
    m_TOldTime = m_T;
    m_TOldNolin = m_T;

    m_Vf.readFromFile(filename_vf);
    m_VfOldTime = m_Vf;
    m_VfOldNolin = m_Vf;

    m_Psi.readFromFile(filename_psi);
    m_PsiOldNolin = m_Psi;

    m_Omega.readFromFile(filename_omega);
    m_OmegaOldTime = m_Omega;
    m_OmegaOldNolin = m_Omega;

    //****************************************************************************************************************************
    //****************************************************************************************************************************
    m_F.resize(m_Net.numOfPoints());

    QVector<int> igl, jgl, igu, ngu;
    fem::triangles::formMatrixProfile(m_Net, igl, jgl, igu, ngu);
    m_A = linal::CSRMatrix(m_Net.numOfPoints(), igl, jgl, igu, ngu);
}
